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Abstract 

We model the dynamics of point Josephson junctions in a ID mi- 
crostrip line using a wave equation with delta distributed sine nonlinear- 
ities. The model is suitable for both low Tc and high Tc systems (0 and 
TT junctions). For a single junction in the line, we found two limiting be- 
haviors: the ohmic mode where the junction acts as a pure resistor which 
stops waves and separates the cavity and the junction mode where the 
wave is homogeneous throughout the strip. This classification allows to 
bound the IV curves of the system. Two junctions in a strip give generally 
ohmic modes and combined junction/ohmic modes and yield information 
about the behavior with an array with many junctions. Finally we use 
this analysis to understand the many junction case for and n junctions 
and the effect of an external magnetic field. 

1 Introduction 

Type I Superconductors are characterized by the phase of a complex number, 
the macroscopic order parameter. Josephson ^ showed that when two such 
superconductors are very close, electrons can tunnel from one to the other ac- 
cording to the relations 



where V and / are respectively the voltage and current across the barrier, (f> is 
the difference between the phases in the two superconductors, s is the contact 
surface, Jc is the critical current density and $o is the flux quantum. The two 
Josephson relations together with Maxwell's equations imply the modulation of 
DC current by an external magnetic field in the static regime (SQUIDs) and 
the conversion of AC current into microwave radiation. 

The ratio of the electromagnetic flux to the flux quantum gives the Josephson 
length, the typical distance of variation of the phase Aj = Depending 
whether their length / is larger or not than Aj, Josephson junctions exhibit 
significantly different behaviors. Small junctions for which I < Xj are used as 
SQUIDS to measure magnetic fields in the static regime or as signal mixers 
where the steep IV characteristic at the gap is used to downshift an unknown 
high frequency signal by combining it with a known local oscillator signal. In 
this latter context Salez et al U suggested to use several unequally spaced 
junctions , instead of just one, to get a better frequency response over a large 
frequency band. Long Josephson junctions for which I » Xj are described by 
the ID sinc-Gordon equation 

0tt - <t)xx + sin(0) = 0, (2) 

which admit the kink solutions varying from to 27r. Long junctions have a 
characteristic dependence of the maximum critical current in the presence of 
magnetic field because of the screening of the applied field by the supercurrent. 
In the dynamical regime long junctions are used as microwave generators but 
provide small output power and poor impedance matching to the transmission 
line. To improve on this, series arrays of junctions have been used, with success 
since they have been shown to synchronize in the presence of external forcing 
[T5] [T^ giving rise to an N'^ power output for N junctions. Such behavior has 
been understood both from the quantum point of view ^2 ^.nd the Resistive 
Shunted Junction (RSJ) equivalent circuit model The cavity surrounding 
the junctions stores radiation which can increase the output power of the device. 

The mathematical description of such a system is a wave equation for the 
cavity coupled to the 2D version of the sine-Gordon equation (0) via interface 
conditions. In one of the authors studied this coupling for a long junction 
with a small surrounding passive region. When the cavity is present only at 
the junction tip and the miss-match at the interface is not too large, the kink 
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crosses it unaffected and adapts its speed. When the passive region surrounds 
the junction, it sets the kink speed. This baUistic kink motion wiU be absent 
for small junctions so in the present work we consider the coupling of a large 
cavity to one or a few small junctions. We want to understand how junctions 
can talk to each other via the surrounding passive region. Two approximations 
are assumed, we consider the ID limit where transverse modes of the microstrip 
do not play a role and neglect the phase variation inside the junctions, in other 
words we neglect the junction length. This is valid for an array of short junctions 
whose length remains smaller than the wave length of the radiation. 

This approach can be applied to high Tc superconductors where due to a dif- 
ferent symmetry of the wave function of the Cooper pairs one has current phase 
relations that involve the higher harmonics [7] . This different symmetry can also 
induce a negative coupling in ^ leading to so-called n junctions. Normal (0) 
and TT junctions can be associated to form O-tt junctions which have semifluxon 
solutions that vary between and tt and are attached to the discontinuity. In 
the static regime the maximum DC current going through the system is strongly 
modulated for a short junction like for a SQUID while it is constant for small 
field when the junction is longer The IV characteristics give resonances 
that are similar to zero field steps W, ^ . Note however the difference between 
the ballistic motion of a kink in a (0) long Josephson junction and the hopping 
of semi-fluxons in a (O-tt) junction. As for junctions, (O-tt) junctions can be 
organized as arrays |1(JII1H and these have interesting paramagnetic behavior. 

The model we introduce allows to describe accurately the coupling between 
small junctions in such an array. Both and tt junctions can be considered and 
disorder can be added. For simplicity we have considered the large damping 
limit where it is possible to develop the analysis. We studied in detail the cases 
of one and two junctions that remain fairly simple and used our conclusions 
for the case of many junctions. For a single junction, we found two limiting 
behaviors: the ohmic mode where the junction acts as a pure resistor which 
stops waves and separates the cavity and the junction mode where the wave is 
homogeneous throughout the microstrip. After introducing the model is section 
2, we consider these limiting behaviors in sections 3,4 and 5. These allow to 
understand in detail the features of the IV characteristics (section 6). The 
behavior of two junctions shown in section 7 follows generally ohmic modes 
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and combined junction/ohmic modes. Using this, in section 8 we generahze the 
study to many junctions embedded in a microstrip and understand the influence 
of the sign of the current density (0 or tt junctions) and external magnetic field. 
Conclusions are presented in section 9. 

2 The model 

The device we model shown in the top panels of Fig. ^ (top view on the left and 
section on the right) is a narrow microstrip in which are embedded a number 
of small junctions in the window design. We proceed from the description 
of the two superconductors as two inductance arrays connected by capacitive 
elements (resp. and a resistor and nonlinear Josephson element) in the passive 
(resp. junction) region. The equivalent circuit is shown in the bottom panel of 
Fig. n In the following subsections we show that the continuum limit yields 
an inhomogeneous 2D sine-Gordon equation which we reduce to ID when the 
microstrip is narrow. For small junctions we then justify our delta function 
approach (^0)) . 

2.1 RSJ model for distributed Josephson junctions 

A simple mathematical model of the window Josephson junction is to describe 
each superconductor by an array of inductances L. The coupling elements be- 
tween two adjacent nodes in each array are, a capacitor C, resistance R and 
Josephson current Ic [2]. The Kirchoff laws at each couple of nodes {ib,it) in 
the bottom and top superconducting layers can be combined to give the relation 
expressing the conservation of currents at node i in the device 

a*..i:^^.A-»|j-|:^». (3) 

where $ = $t — <!>(, (resp. = — '^b), is the phase difference between the 
two superconductors in the junction (resp. passive) part, the summation 
is applied to the nearest neighbors, /f is the critical current which is non zero 
only in the junctions and Ri is the resistance associated to the current of quasi- 
particles again non zero only in the junctions. In agreement with experiments, 
we have assumed the same surface inductance for the whole sample but have 



4 



taken into account the variation of the capacity C and resistance R from the 
junctions to the rest of the cavity. Note that equation Q is a discretisation of 
Maxwell's equations (wave equation part) and Josephson constitutive equations 
(sine term), assuming an electric field normal to the plates, a magnetic field in 
the junction plane and perfect symmetry between the top and bottom super- 
conducting layers. We can now obtain the model in the continuum limit, more 
suitable for analysis. 



2.2 Continuum limit 



The continuum version of the system |(2J) can be derived by introducing the 
following quantities per unit area {a?) of elementary cells of length a. 



C. 



J 



9 ' Ji o ' 



We normalize the phases by the flux quantum $0, 

$0 

and introduce the Josephson characteristic length 



A5 = 



(4) 



(5) 



(6) 



Notice that in this 2D problem the inductance associated to each cell is equal to 
the branch inductance L. To differentiate the junction regions from the passive 
regions we introduce the indicator function gi such that gi = 1 in the junctions 
and elsewhere. We can then write 



9if , =^giz: and C^ = giCj + {1 - gi)Ci 
n r 



(7) 



where j'^, r, Cj and Ci are respectively the critical current densities, the quasi 
particle resistance in the junction, the capacitance per unit area in the junction 
and the capacitance per unit area in the passive region. 
We substitute relations Q - I7J in and obtain 



L{Cj - Ci)(pi + simp., 
A -, 



r 



= 



which in the continuum limit yields 



LCiiptt - A(p - 



LX^j{Cj - C/)(ptt + sin ^ 



= 0. 



(8) 



(9) 
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To obtain the final equation we rescale time by the inverse of the plasma fre- 
quency Loy^ = Xj\/ LC I so that i — ujjt. As unit of space we use the Josephson 
length, X = x/Xj, y = y/Xj. The final equation is then 

(fitt - Alp + g{x,y){c(ftt +sm(p + aipt) = , (10) 

where the tildes have been omitted for simplicity and where the coefficients c 
and a are 

— — -1 — (II) 

^~ Ci ' r ^Ci ^ ' 

The boundary conditions represent the input of an external current / or 
magnetic field b on the device. They are 

(12) 

where < < 1. The case v = \ corresponds to a pure inline current feed and 
= to a pure overlap current feed. We will mostly consider the latter case in 
the discussion. 



2.3 Reduction to a ID problem 

When the microstrip is narrow and the junctions occupy the whole width we 
can assume that g is independent of y. We can then write 

2 ^ 

<^(x, y, t) = + 0„(x, t) cos (^) , (13) 

n— 

where the first term takes care of the current feed and the second satisfies 
homogeneous Neumann boundary conditions for y — 0, w. 

The calculations detailed in |51 show that for a narrow width and small 
current only the first mode (p^ needs to be taken into account. We then obtain 
the following equation for <j)Q where the has been dropped for simplicity 

<Att - 4>xx + g{x){c(f)tt + sin(/) + acf)t) v- , (14) 

s 

with the boundary conditions 4'xx=o = ^(1 ^ 't>xx=i — ~ ^)ts- ^^"^ 

following we will write 

^■ = ? (^^) 
and assume an overlap current feed (i^ = 1) except in section 7. 
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2.4 The delta function model 

As the area of the junction is reduced the total super-current is reduced and 
tends to zero. To describe situations where the area is small so that the vari- 
ations of the phase can be neglected in the junction but the supercurrent is 
significant we introduce the following delta function model. Introduce the func- 
tion g{x) = dagh{x) 

l/{2h) a~h<x<a + h 
elsewhere 

(t^tt - (l>xx + dagh{x) (sin(0) -I- a(f>t + apu) = j- (16) 



9h{x) 



In the following we assume to be in x and t and (jittx to be continuous. 
We now show that equation Ijltill converges when — > to 

<^tt - <l)xx + daS{x - a) (sin(0) + a(j)t -f c(j)tt) = j- (17) 

For that introduce the Taylor expansion of (f) near x = a 

(j){x, t) = cj){a, t) + {x - a)(j>x{a, t) + 0{{x - af). (18) 

In the following for simplicity of writing we will omit the t and x dependencies. 
We integrate from a — h to a -\- h (see details in appendix: Continuum 
limit). The first term 

4>ttdx — !■ O when /i — > 0, 

-h 



the second one 



and the third one 



a+h 

(j)xxdx ^ cf)x{a + h)- (j)x{a), 

-h 



da 



h^o 2h 



a-\-h 



lim — / (sin(0) + a4>t + c4>tt)dx = da (sin(0(a)) -f- a(j)t{a) + ccjjttia)) , 



-h 



when /i — !■ which is consistent with the model 

(t)tt - <f>xx + daS{x - a) {sm{(t>) + a(j>t + ccj)tt) = j, (19) 
together with the boundary conditions (a; = 0) = , (px {x = l^t) ~ Q. 
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In the following we will consider the case c = for simplicity. We assume 
the same capacity in the junction and in the passive region so that the final 
equation is 

4'u - (f>xx + da6{x - a) (sin(0) + a(f>t) = j. (20) 

This equation can be considered as a ID model of a point Josephson junction 
embedded in a microstrip cavity. 

3 The single junction case 

The partial differential equation H2(J|) involves a distribution, however its solu- 
tion is differentiable except at the junction points. To solve it numerically we 
integrate the equation over reference intervals using the finite volume approach. 
The resulting system of ordinary differential equations is then solved using the 
D0PRI5 4th and 5th order Runge-Kutta variable step integrator developed by 
Hairer and Norsett ^7] . Details of the numerical code are given in the second 
part of the appendix. 

We now consider the two limiting cases of a static solution and high voltage 
for which simple expressions of the solution can be obtained. In particular these 
have been used to validate our numerical procedure. 

3.1 Static behavior 

We assume 0t — (pu = in and get 

- (I)xx + S{x - a)dasm{(f)) = j (21) 
'px\x=o.i=0 (22) 

For a; 7^ a, —(jyxx = j so that for x < a, <!) = 4>i = —^x'^ + Cix + Di and for 
X > a, (j) = (j)r = ~2^^ ^'■^ Integrating H21|l over the whole domain we 

get the conservation of current 

dasin{^ia))^jl, (23) 

so that no solution exists for \jl\ > da- 

The solution is obtained as usual by assuming continuity of the phase and 
using the jump condition on the gradient obtained by integrating on a small 
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domain including the junction 




da sin(0(a)). 



(24) 



Using these two constraints we obtain the static sohition 



4's{x) 



i(a — x){a + x) + arcsin i^j^j 
|((a x){x + a + 21) + arcsin 




X < a 



a < X 



(25) 



3.2 High voltage behavior 

Since the Josephson current is bounded, its influence in (|2U|) will decrease as the 
current j is increased. In this limit of high current the phase is growing very fast 
so that (pt is close to the average voltage V = {(t)t)f We can write approximately 
(j){x,t) — Vt + (j)v{x) so that we can average (pn|l assuming {(j)t{x,t))^ — V, 
(sin((/)(a)))j = and {(j)tt)t = 0, yielding 



with the boundary conditions (t>vx\x=o.i = 0- Integrating this equation over the 
whole domain we get V = -f—. 

Using the same constraints as in the static case (continuity of the phase and 
jump of its gradient at the junction) we get up to a constant the high voltage 
solution (j)^ 



Notice that this and precedent solutions have the same spatial behavior; (ps — 
4>y = Vt + C. We will not find this result with more junctions (except in 
particulary cases). 

3.3 Fourier representation 

Because of the homogeneous Neumann boundary conditions one can write the 
solution as a cosine Fourier series 



4>yXx + 6 {x — a) daOiV = j, 



(26) 




(27) 




(28) 
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The solution of (|20|l has a discontinuous first spatial derivative to that the 
amplitude of the mode An decreases as 1 /n^ jJS] . This means that in general 
many modes are needed for the description of the solution. However since the 
system is close to a linear one we will see that these modes provide insight into 
the limiting behaviors. 

Plugging into (pn|l and projecting we get the evolution of the modes, 

lAott + da{sm{(l)a) + a(j)at) - jl = 0, (29) 

— j A„ + c„ (sin(</)a ) + a0at) = 0, (30) 

where we have introduced the coupling coefhcients associated to the junction at 
X = a, Cn — cos (^^). Equations and l|5n|) are coupled through the term 
sin((/)a) + a(l)af Due to this particular feature, two limiting behaviors can be 
expected. 

• The ohmic mode for which 

daC(4'a\t - jl = . (31) 

Then the equations for the An have a right hand side so that the cavity 
modes are driven by the junction. 

• The junction mode such that 

daism{(t)a) + a(l)at) ' jl = ■ (32) 

In this case the equations for the An have a constant right hand side 
{— jl). We will show that in this case the cavity is synchronized with the 
junction. 

The equations giving the evolution of (/)a in these two limiting behaviors can 
be solved yielding for the ohmic mode 

(ka = -y-t + Ci , (33) 

"a a 

where Ci is a constant and a voltage 



K = /-. (34) 
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For the junction mode the equation can be solved using separation of variables 
and we get 



2 arctan i 



where the junction mode voltage Vj is 



aFjtanl ^(t + Cs) ) +1 




(35) 



(36) 

ct y \ "a / 

where C2 is a constant. 

Fig. |5]shows these two limiting behaviors for a junction placed in a = —1.57 
with da = 1 Z = 10 and a current j ~ 0.1374. The numerical solution is 
presented with the crosses while the analytical estimates (|33|l and H35|) are shown 
in continuous line. The left panel corresponds to the ohmic mode and one can 
see the excellent agreement of the numerical solution with the line given by (|33() . 
The right panel presents the junction mode and the extension of the estimate 
(see below), again the agreement is very good. Notice that for the same time 
interval the variation of 0t(a,t) is different for the two panels of Fig. [3 so that 
we will get a different average voltage for these two behaviors. 



4 The ohmic mode 

In this case (/)t{a,t) = Vo = To get the field we need to solve on each 

subdomain [0, a[ and ]a, I] the following problems 

{<l>tt - 4>xx = j J 4>tt - 4>xx = j 

Consider the left subdomain [0, a[, we introduce (p {x, t) = v (x, t) + so that 
if 4> solves the left problem then v solves 

vtt - Vxx = j (37) 
^^.lo = vtl = (38) 

A particular solution of H37|) is vi = — |a;^ + Di where Di is an integration 
constant. Then w = v — vi solves the homogeneous wave equation 

dttw - dxxW = 0, 



11 



with Wx\o = and Wt\a — 0. We can assume without loss of generaUty that 
w (x, 0) = 0. The condition Wx\q = imposes a solution of the type 

w {x, t) — ai sin (Vji) cos (Vix) , (39) 

where Vi is a constant. Also, wt\a = imposes a condition on Vi because 

Wt {x — a,t) — ai sin (Vit) cos(V/a) — 

is possible for all i's only if cos (VJa) = 0. So that we get 

2kl + lTT 

Vl = — ^ ■ 

2 a 

The problem is identical for x £ ]a, I] and using similar arguments we get Vr = 
i 

2 



[x — I) + Dr and 



2kr + l TT 

Vr — 



2 a- I 

We obtain a solution for the equation (|20|l . The phase on the whole domain is 
then 



(l>{x,t) 



+ ai sin (Vit) cos (Vix) + ^{a — x){a + x) x < a 

+ ar sin (Vrt) cos (V^ (a; - 0) + ~ ^)('^ + a; - 2^) a < a; 

(40) 

From the problem (|20|) integrating the equation on a small interval centered on 
X = a and using H31|) we get the jump condition 

+ jit 
[4>x]a- - jl = da sin(0(a, t)) = da sin(— — ). 

ada 

On the other hand from H4()|l we get 

[(j}x]a- ~ jl = ai sin (Vit) sin (Via) — GrVr sin {Vrt) sm{Vr (a — I)) 
= Ci sin (Vit) - C2 sin (Vrt) 

Since Ci and C2 are independent of t, this implies that sin{Vit), sin{Vrt) and 
sin(^j-) have the same period. We obtain, 

2fc7 + 1 TT jl 2fcr + 1 TT , ^ 

— —- = = — — = V (AI) 

2 a ada 2 a~l ' ^ ' 

if ai and are not equal to zero. Then one can write the solution as 

j 4>y{x,t) + ai sin {Vt) cos {Vx) x<a 
(p[x,t) = < (42) 
I 4>y{x,t) + ar sin {Vt) cos {V {x — I)) a < x. 
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left only: (j){x, t) = 
right only: t) 



where is the high voltage solution. Because of this the instantaneous voltage 
4)t has a discontinuous first derivative at 2: = a. 

Assuming that a = 1 we see that the consistency equation H41f) implies 

(a — l)2ki — 2akr ~ I. 

There are very few solutions for a rational et I integer and for most a there are 
none. Let us look at other possibilities. Notice first that space homogeneous 
functions (a; = or = 0) are solutions. Having said this it is easy to see that 
one can obtain left only or right only oscillations. 

(f)y{x,t) + ais\Q-{Vt)cos{Vx) x<a 
(f)v{x, t) a < X 

I (/)t,(x, t) X < a 

[ (j}vix,t) + arSin{Vt)cos{V (x — I)) a<x 
Again the instantaneous voltage at the junction has a discontinuous first deriva- 
tive. 

Fig. 13 shows the instantaneous voltage as a function of x for different times 
for two such ohmic modes, the left only oscillation (left panel) and right and 
left oscillations (right panel). The left panel shows (ptix, t) for 60 values of time 
equidistributed between and 10. The parameters are the same as for Fig. El 
In this case we have fc; = 1 so that about 1.5 periods exist to the left of a; = a. 
Notice how the amplitude of 0* is practically constant and equal to the average 
voltage V = jl = 1.37 on the right of x = a. The right panel corresponds 
to a right and left oscillation for a = 2 and we show three consecutives times, 
t ^ 0, 0.125 and 0.250. Here = so that we get about 1 period to 

the left of a; = a and two periods to the right of a; = a. We show for comparison 
V = ^ {— jl) the average voltage. Notice in both cases the discontinuity of (pt 
at X = a. 

Let us show that this is consistent with the whole equation (|20|l . for the 
intensity I — V/ ada ■ Consider the left only solution 

r Vt+f sin(yt) cos{Vx) - I [a;2 - a2] x<a 
^ \ Vt~^[{x-l)^ + ia~l)^] x>a 

with, V — ^^^^f • On the interval [0, a[ (respectively ]a, I]) (j) satisfies the linear 
wave equation (jjtt — 4>xx — j with the boundary condition <f)x\o — (respectively 
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= 0). We now integrate (|2U|1 on [0, /] to obtain 

pa pi pi 

I (j)udx + / (t)udx + / (/Jxxdx + da sin (Vt) + adaV = jl, (44) 

Jo Ja Jo 

SO that 

-ai sm{Vt) [sin(Fa;)]o + + + rfa sin = 0. 
The end resuh is 

{-if'+'^ai siniyt) = da sm{Vt) 
which impUes a; = (— 1 ^ da and the equation is balanced. 

5 The junction mode 

The expression (|35f) a priori defined only for a given time interval can be C°° 
extended for all times using the continuation 

where k is an integer. The pseudo-period T oi (f)a is T — , ^"^ so the 

voltage isV—^ — j^ \/ {j/daY ~ 1- This can also be written 

J = y\/(«^)' + l (45) 

So we deduce from A^^ (t) = and the voltage expression ^ that Ao{t) = 
Vt + C. 

Let us now calculate the contribution of the Fourier modes and give a solution 
to the whole equation l|3U|l . From H3U|) and the previous equality, we obtain: 

^;W + (^)'^pW = -2jCp (46) 

The solution of this equation, Ap is: 

Ap (t) = -fp cos (^<) + f]p sin (^t) - ^^cp (47) 

Now we have the expression of all the Fourier modes. We know that 0(a, t) = 
Y^n=o ^n(^)'^n- Substituting the Fourier modes (|47|l into this expression we 
obtain: 

-t-oo 

0(a,i) ^Vt + C + J2 
p=i 
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The series in the last term is uniformly convergent so that the term can be 
collected in the form of a constant yielding the final simpler expression 



+00 



(a, t) = Vt + C2+Y^ [ipCp cos (^i) + ppCp sin (^t) ] (49) 



p=i 



Introduce (/)(a, t) — Vt + C2 ^ u{t). We can choose C2 as /J" u{t)dt = 0. Notice 
that It is a C°° and T-periodic function. We can make Fourier projection on it. 



Until Cp 7^ 0, we have 



i(i) = ^ \lpCp cos {^ij + PpCp sin {^t 



2 '■^ 



"fp = ^]^ cos i^t) dt (50) 

Pp = -^ I <i) sin i^t) dt (51) 
Cpi Jq \ I y 

We can make the following remarks 

1. V =^ ^ ^ so T ^ Yi and w is |^ periodic. 

2. If = ^ the terms of the Fourier series that aren't zero are the multiples 
of n. 

3. If Cn = then the junction modes for V — ^ can't exist. For example 
when the junction is in the middle of the circuit, then c° = for n = 2fc + l 
and junction modes exist only for V — 2^ where k is an integer (see left 
panel of figure [SJl . 

4. When \cn\ is very small, the circuit in the junction mode stores a lot of 
energy but is unstable. We will see this in the next section. 

Now remember (j28|l . we can give the explicit solution of the phase in the n th 
junction mode, up to a constant 

+00 

, , , mr v-^ r /pmr \ . /pmr \l /pmr \ 

(t)(x,t) ^ —t + ^^'jp cosy-j— tj + Pp sm\^— l—tj ^ cosy— j—xj . (52) 

p=i 

Notice that the instantaneous voltage (f>t has a continuous derivative at the 
junction position x = a. Fig. ^ shows the instantaneous voltage (j)tix,t) as a 
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function of x for 60 consecutive times for the same parameters as in Fig. |21 
Notice how the solution oscillates over the whole domain and that one cannot 
see the junction position. This type of motion is similar to the transparency 
observed in the reflection coefficient of a stripe ^Sj when an integer number of 
half-periods of the wave " fits" in the stripe. 



6 Single junction case: IV characteristics 

For a large enough time the system finds a stationary state where the energy 
provided by the input current is balanced by dissipation. This state is described 
by the IV characteristic curve. This is measured on real devices and it is there- 
fore very important to understand its features. 

To compute it we fix the current and run the system for about ti = 3000 
time units so that all transients have died off. Then we compute the average 
voltage 

/ , / ,^^ _ (p{a,ti+T) - (j){a,ti) 

{Ma, t))t = Tf, (53) 

at the junction. The value T is about 2000. The space average of the average 
voltage in time (^j (j)t{x,t)dx^ is compared to this quantity to ensure that 
the system is completely thermalized. Only then is the voltage recorded. 

The IV curves are computed for 99 steps of increasing current j from to 
4 starting with an initial condition (f){x,0) = and ^((a::, 0) = 0. We then 
decrease the current from 4 to 0. All plots show these two curves for increasing 
and decreasing current. 

6.1 Bounds for IV curves 

We have isolated two characteristic behaviors of the circuit. These allow us to 
bound the IV curves for a single junction in a microstrip. Fig. [S] shows two IV 
curves for a centered junction a = 5 (left panel) and an off-centered junction 
a = 9 in a microstrip of total length I = 10. For each case we plot the limiting 
behaviors given by the ohmic mode and junction mode voltages ()34|l and (|36fl . 
V, j are bounded by the relations 



daaV <jl < da^iaV)^ + 1 
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The ohmic modes, (resp. the junction modes) only occur for voltages such 
that V ^Vi = ^V^f or F = K = (^esp- V = ^). The vohage 

(horizontal) axis is labeled with the indices corresponding in each case to the 
junction modes. The ohmic modes are marked with ticks but not labelled. 
We always go from an ohmic mode to a junction mode and vice versa. As the 
current is increased the junction mode ceases to exist and the system then jumps 
to the closest ohmic mode. Notice the hysteresis obtained between the curve 
for increasing j and the one for decreasing j. Both the junction voltage and 
ohmic mode voltages arc fixed by the length of the system. On the other hand 
the current for which one obtains these voltages depends on the damping a and 
strength of the junction da- For fixed da, an decrease of a will tilt the IV curves 
towards the horizontal voltage axis and this will separate the cavity modes and 
cause large voltage jumps from a cavity mode to a far away ohmic mode and 
hysteresis. Conversely an increase of a will make the IV curves more vertical so 
that the voltage jumps and hysteresis will be reduced. Note also that as in the 
static case, a magnetic field has no effect on the IV curves for a single junction. 

When the junction is centered as in the left panel of Fig. Othc ohmic modes 
(resp. junction modes) correspond to the odd V — {2k + l)7r// (resp. even 
V = 2kn/l) cavity modes so the voltage interval between a junction mode and 
its corresponding ohmic mode neighbors is constant. This does not happen 
when the junction is off-centered as in the right panel of Fig. |S1 so that one 
can get very sharp resonances connecting an ohmic mode to its junction mode 
neighbor like the ones for n — 2,S and 4. The resonance for n = 3 was obtained 
by decreasing current. For larger values of n the voltage separation becomes 
larger so that the resonances are softer. 

6.2 Study of resonances 

We have been able to bound the IV curves for our system. Now let us discuss 
the fine structure of the resonances, why some are sharp (see resonances 3,4 and 
6 in Fig. [SJ while others are not. To explain this one should consider the energy 
of the system and the Fourier modes. Multiplying H2U|I by 0f and integrating 
from to I we obtain the work equation 



if 

2 Jo 



+ (jj^dx + da (1 — cos (/)(a, t)) 



j (t)tdx - daa(t)^{a,t). (54) 
Jo 
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We can write ^ as dfEt = dt{Ep + Ej) = Pdc + Pq where Ep ^ \ J^{(j)f + 
(j)^)dx is the energy in the passive region, Ej = (1 — cos (j){a, t)) the Josephson 
energy, P^c — j jl^(j)tdx is the power given by the direct current and P^ = 
daa(j)f{a,t) is power dissipated by the quasiparticles. In figure we plot from 
top to bottom Et-, Ep, Pq and Ej for a centered junction a = 5 (bottom panels) 
and an off-centered junction a — 9 (top panels). The left panels show the 
junction mode n — A and the right panels the corresponding ohmic mode for 
the same current. In the junction mode the quasi-particlc term Pq oscillates 
strongly and the Josephson energy Ej is very anharmonic. In the ohmic mode 
Pq is almost constant and Ej is sinusoidal. In all cases the energy in the cavity 
Ep is much larger than the junction energy Ej . When the junction is off-centered 
(top left panel) the energy stored in the cavity is much larger. This is due to 
the smallness of the coupling coefficient c„ as we now show. 

For that consider the phase 4>{a,t) at the junction as given by (|35|l . It is 
independent of the junction position a. Since we are in the junction mode 
V — mr/l we write from H52(l 



The first statement implies that the terms jpCnp and (3pCnp are constants inde- 
pendent of the junction position a so that jp and Pp are proportional to l/c„p. 
As a consequence if c„p is small the junction mode of index n corresponds to a 
large amplitude in the mode Anp. The system accumulates this energy in the 
passive region and the resonance is very sharp. 

Fig. El shows a plot as a function of time, in each panel, from top to bottom 
of the total energy Et, the energy in the passive region Ep, the power dissipated 
by the quasi-particles Pq and the Josephson energy Ej for a centered junction 
a = 5 (bottom panels) and an off-centered junction a = 9 (top panels). The left 
panels show the junction mode n — A and the right panels the corresponding 
ohmic mode for the same current. The coupling coefficient is small C4 = 0.3 
in the top panel while it is maximum C4 = 1 in the bottom panel where the 
junction is centered. Therefore one expects a sharp resonance and a highly 
excited passive cavity for a = 9 as opposed to a = 5. This is indeed the case 
as one can see from the plot of Ep in the left bottom and left top panels. The 
Josephson energy Ej and the power dissipated by the junction are of the same 
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order. In the ohmic regime this accumulation of energy in the cavity is absent. 
There is a small difference in Ep between the two configurations because they 
have different lengths of passive regions. Note that the total energy is practically 
independent of time and so is the dissipation. There is an exact balance between 
the Josephson energy and the energy in the cavity Ep. 

7 The many junction case 

Assuming a general situation with magnetic field, a mixed current feed and 
many junctions in the micro-strip we can generalize our model to 

4>tt - 4'xx + '^ djS{x - aj)(sin((/)) + a(j)t) = vj (55) 

'j 

with the boundary conditions (t>x\x=o — h—'^{l^v) and (t>x\x=i = b+ ^{l — v). 
If = 1, In the model (|55|l ly = 1 (resp. 1^ = 0) corresponds to an overlap (resp. 
inline) current feed. 

7.1 Static case 

We have seen that for a single junction the maximum current jl = da is reached 
for any value of the magnetic field. This is not true for two junctions (a SQUID) 
where the maximum critical current if obtained for 6 = mod. where 

02— ai 

n is an integer |22j . Static solutions for many junctions in a ID microstrip will 
be discussed in a forthcoming paper [231 ■ 

7.2 High voltage case 

As in the case of a single junction, when the stationary state is reached we can 
assume (l)t = V and average out the (f>tt and sin(0) terms yielding the following 
boundary value problem for the time independent part (j)y of 4> 

n 

- (l>vxx + ^dkS{x - ak)aV = vj, (56) 

k=l 

</>.Jo.i =&T^(l-i^), (57) 

On each interval separating the junctions (t>vxx — ^'^ that 0„ is a second 
degree polynomial 4>v{x) = — (i^j/2)x^ + bx -\- d. 
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Integrating (|56|) from to I yields the voltage 

Integrating (|56|) over small domains containing the junctions we obtain as usual 
the jump conditions in the derivatives 

Using these remarks we can calculate (f)^ for any number of junctions (for 
simplicity one can assume 0„(ai) = 0). Let's introduce S{x), the spatial curva- 
ture of {(t>)^. For two junctions we get 

{^{x^-al) + (b-'^{l-iy)'^ix-a,) x < ai 

^{x^ - «?) + +&-#(!- ^)) - «i) a,<x< a2 

^{x^ - al) + (b- f (1 + t/)) {x - ai) + ^2 a2<x 

(58) 

where S2 = ^(ai " a?) + (diq^ 6- ^(1 - t^)) (02-01) is the phase differ- 
ence between the second and first junction (S2 = 4'i'^2, t) — 0(ai, t)). The phase 
is then given by 

(l){x,t) ^(t)^{x) +Vt + C, (59) 

manyj where C is an integration constant. Notice that it is possible to find 
4>v{x) for an arbitrary number of junctions and compute the phase difference 
4>v{ai) = Si = cj){ai,t) - (t){ai^i,t). 

In Fig. [7| we show (/)(x, t) for 6 consecutive values of time separated by 
St = 0.1 for a system of two junctions located respectively at oi = 2 and 
02 — 6.2 in a domain of length / = 10 with current j = 1 fed in overlap 
geometry (i^ = 1). One can see that the behavior is close to what is predicted 
by the above expression. 

7.3 Fourier representation 

We can use the high voltage solution to simplify (|55|l . For that take i/j — 4> + 4>y 
so that H55() becomes 

iptt - i^xx + dk6{x - Qk) sin('(/' - Sk) + aipt - — t ) = 0, (60) 
fit V Ek=idkJ 

Mx^o-j = 0- (61) 
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The transformation using the high vohage solution concentrates the current on 
the junctions. In this representation inhne or overlap configurations differ by 
the terms Sk- 

Let us go back with Fourier series. Because of the homogeneous Neumann 
boundary conditions on Ht)U|) . we can decompose -0 on a cosine Fourier series 

+00 

n=0 

and obtain the following equations for the modes An 

Mo + V fsin(7/.j -Sj) + a^jt - , ) = 0, (62) 

where c:^ = cos ( ) , -ifij = ^jj{aj,t), ipjt = dti^{aj, t) 



8 IV curves for two or more junctions 

8.1 2 symmetrically placed junctions, limiting behaviors. 

In a symmetric circuit, we have 02 = I — ai and di = d2- For simplicity we 
assume no magnetic field so that 6 = too. With this symmetry, it is easy to 
show that C2J. = Cjj. and c^^+i ~ ~'^2fc+i- important to remark that when 
6 = 0, 5*2 = 0. So there is no phase shift between the two junctions. Let us now 
go back to the Fourier modes, H62() becomes: 

di (sin(V'i) + aV'it + sin(V'2) + aV'2t) - = -^^o- (64) 

For p — 2k, 

II /2fc7r\^ 2c} 
^fe+( — ) = p[di(sin(V'i) + aV'it +sin(V'2) + ai/'2t) -jZ] 

We want to show that the following equation for A2k+i is decoupled from the 
system of Fourier modes 

^(2fc+l)^V^ _ -2cife^i 



-42fc+i + ^ j j ^2fc+i = j di(sin(i/ii) - sin(V'2) + a{^it - '02t))- 

(65) 
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For that we group together the even terms and odd terms of the Fourier series of 
tpit a,ndip2t and define So (t) = J2t=o^2k^2k{t) and Se{t) = J2k=o^2k+i^2k+i{t)- 
Notice that: 

•0(01, i) = So(i) + Se(t) and ^{a2,t) ^ Soit) - Se{t) 

Replacing ipi and V'2 in equation (|65|l . we obtain: 

„ ^ (2fc + l)^ y ^ -4ck.+i , f f , ■ ( ^ '\ laa, 
^2fc+i + I ^ I ^2fc+i = ^ di (^cos(se)sm(so) - asgj (66) 

This equations show that the system performs like a single centered junction 
with different coupling coefficients. So we can balance the Fourier modes in the 
junction mode: with initial conditions V/c, A2fe+i(0) = 0, ^2fc+i(0) = we have 
5*2(0) = and as a consequence A2k+i{t) = 0, G 3?+. This shows that the 
equations for the odd terms are decoupled from the system. 

Junctions modes. 
As we see, the system (|^ behaves like a single junction: the odd terms are 
equal to zero and the even terms give the junction modes (see H35(l ): 

fdi r /I 

01 = 09 = 2 arctan < — 



tan( — y,(t + C2))y, + l 



Notice that the two junctions are perfectly synchronized. In a non-symmetric 
situation, junctions modes for both jimctions are particular cases. 
Ohmic modes. 

Choosing 6*2 = allow us to show again analytical solutions. Assuming 6 = 0, 
there are three types of solution, given by different oscillations of the phase: 

1. The phase oscillates between the junctions so that V — -^l^r^^. Due to 
the symmetry of the problem, we cannot find solutions with V ~ al'^ai ■ 

2. The phase oscillates outside the junctions: V — -^^^^-p^- 

3. The phase oscillates in all the passive parts so that V = ^^^^^^^^ = 
This is only possible if 02 — ai is a fraction of / — 02. 

To conclude in the symmetric case, we know limit IV curves with exact solutions. 
See the left panel of Fig^] In the general case, we cannot find exact solutions, 
we have near ohmic mode solutions and peaks at the expected resonant voltages 



22 



8.2 The non-symmetric cases: IV curves 

The IV curve obtained for a circuit of two junctions placed in a non symmetric 
fashion confirms the existence of ohmic modes at the expected voltages. It is still 
possible to characterize the IV curves using ohmic modes and junction modes 
even if there are only approximate. Without symmetry and when 5*2 = [tt] 
we can obtain exactly any one of ohmic modes discussed in previous section. 
In addition to the symmetric case there two new possibilities where V can be 
computed: 

1. The phase oscillates between the junctions and to the right of ai. See the 
left panel of Fig|Slas an example. 

2. The phase oscillates to the left of a2 and between the junctions. 

In Fig|Sl the middle panel shows the second point of the symmetric case and the 
right panel illustrates the case where the phase oscillates between the junctions 
for the non-symmetric case. 

Peaks appear at the junction modes V = but they are in general smaller 
than in the symmetric cases. Fig El shows the IV curves for two non symmetri- 
cally placed junctions and indeed one can see that the height of the resonances 
is lower than one for two decoupled junctions (the upper curve in the three 
panels). By varying the junction critical current (the coefficients di and d2 ) we 
show that the resonances are linked to one of the two junctions. To understand 
this phenomenon, go back to the coupling coefficients c:^. For the second reso- 
nance in Fig. IC2I = 1, Ic^l = so that the junctions are strongly coupled 
to the system. This is like in the symmetric case but it is an exception. The IV 
curve shows that they are together in junction mode (or very close). Consider 
now the third resonance. The coupling coefficients are Cg = and c\ = 0.90 
so that the junction ai is decoupled from the system, we can assume ohmic 
behavior for it and 02 is in junction mode. The fourth resonance corresponds 
to the reverse situation, c\ — 1, c\ = 0.06 where a2 is decoupled. 

To confirm this, we plot in Fig. 1101 4>i and 02 for the resonances 3 and 4 with 
di — 0.7, c?2 = 0.3 (left panel of Fig. |5J and di ~ d2 ^ 0.5 (middle panel of Fig. 

Together with the numerical results we plot the junction mode solution H35|l 
in dashed line. For the top panel (third resonance) one can see that the junction 
02 follows closely this behavior as opposed to ai which is in ohmic mode. The 



23 



situation is reversed for the fourth resonance as shown in the bottom panels. 
For resonances 2, 3 and 4, the value of the critical current di ,i — 1,2 does 
not modify appreciably the behavior of the phase at the junctions (remark that 
5*2 9^ but stay small). 

For a general resonance the junctions will work in different modes, one will 
be closer to a junction mode while the other will be closer to an ohmic mode. 
We call this mode a resistor-junction mode. 

Resistor-junction modes. 
Resistor-junction modes exist if the circuit has one junction heavily coupled 
while the other is decoupled. In this case, the system behaves almost like a 
single junction circuit. This is why these modes are associed to voltages V = 

If the junction is in ohmic mode, then = dkSm{Vt + Ci) + Ck- 

So that we can only have only one harmonic in the circuit. Whereas a junc- 
tion mode exists only with an infinity of harmonics. So the system reaches a 
constructive equilibrium. For example in all the panels of Fig. 1101 we plot the 
exact junction mode to compare. The observation of the phases in the junction 
mode indicates that it lingers more around (2fc+l)7r/2 than the expression (|35(l . 
This gives for the IV curves of Fig. IHl resonances that are slightly higher than 
expected. This summarizes the influence of one junction on the other. 

8.3 7r-junctions and magnetic field effects 

We now consider the situation where the critical current di for some i can be 
negative. Then the ith junction is a so-called tt junction for which the Josephson 
relation is given by sin(0 -|- tt) instead of sin((/)). In this case we will show 
that the dynamics is very different than the one where all the di are positive. 
Therefore combining tt and standard junctions provides a way to change the 
resonant frequencies of the array. Another important parameter which we have 
not varied up to now is the magnetic field h. We will see that a change of h 
changes the IV characteristic in a way that is similar to changing the critical 
current density of some of the junctions in the array. This is another way to 
tune the array to resonate on specific frequencies. 

We now illustrate these specific points on some examples. Consider an array 
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of N point junctions described by 

N 

4>tt - (l)xx + '^ djS{x - aj){sm{(j>) + a(j)t + c(j)tt) = vj , 

with the boundary conditions (j)x\x=o — b+^{l~v) and (j)x\x=i —b—'^{l — i'). 
First note that if < for all i's then we can take -0 = + tt to get us back to 
the standard case di > 0. We therefore obtain the same IV curve as when all 
the di > 0. 

An interesting case corresponds to — 2. Let us consider the effect of 
an additional magnetic field b' on the IV curve. For that introduce ip{x,t) = 
(j){x,t) + f{x) where f{x) — b'{x — ai). The equation verified by '0 is 

fAtt - i^xx + di6{x- ai) (sin(0 + /i) + ai/^f ) + ^2(5 {x - 02) (sin(?/; + /2) + aipt) = vj, 

with the boundary conditions il!x\x=n = b + b' + ^{\ — v) and '4>x\x=i = b + b' — 
— v). When fi—0 and /2 = 2tt ie when b' = 27r/(a2 — ai) we obtain the 
same equation for and "0 apart from the boundary conditions. The IV curves 
for the two magnetic fields 6 = and b = 2mr/{a2 ~ a-i) where n is an integer 
are then identical. Now if one chooses a magnetic field b — n/{a2 ~ ai) then 
the phases are shifted by tt one from the other and this is like changing the sign 
of one of the di . In practice these results only makes sense for magnetic fields 
for which the size of the junctions can be neglected. This is the limit of our 
approach. 

Fig. ^Jshows two IV curves for two symmetrically placed junctions ai = 2 et 
02 = 8 such that b ^ di — d2 — 0.5 corresponding to standard junctions (left 
panel) and b — di = — ^2 = 0.5 corresponding to a standard junction and a tt 
junction (right panel). We observe in the left panel resonances (junction modes) 
every exactly as for a centered junction. In the right panel, resonances 
(junction modes) are found every (^^LtlliL^ xhc tt phase shift between the two 
junctions implies that odd harmonics are highly coupled and even harmonics 
are decoupled from the system. We could have obtained the same change in the 
nature of the resonances by using the standard array (left panel) and changing 
the magnetic field to 6 = 7r/(a2 — ai). As the field is changed one then goes 
continuously from an IV curve with only even resonances to an IV curve with 
only odd resonances. 
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Now let us consider how these results obtained for TV = 2 affect the behavior 
of an array of A'' > 3 junctions. If there is a d such that o^+i — = rijd where 
Ui is an integer for all i then the IV curves for 6 = will be the same as for 
b = 271/(1. If one chooses b = ir/d the junctions such that fi = n mod. 27r will 
be reversed while the ones for which /,; = 2n mod. 2t: will not be affected. This 
is similar to a flute which is tuned by changing nodes. 

9 Conclusion 

We have introduced a simple and general model using a wave equation with 
delta distributed sine nonlinearities to describe the dynamics of point Josephson 
junctions in a ID microstrip. This can be extended to tt junctions for which the 
ciirrcnt density can be either positive or negative in the domain. We can also 
apply it to situations where the current density in each junction is different. 

For a single junction we have shown two limiting behaviors, the ohmic mode 
where the junction behaves as a resistor driven by the cavity, separating waves 
and the junction mode where the Josephson element is driving the cavity. These 
limiting behaviors have allowed us to bound the IV curves and understand the 
observed resonances. When another junction is present in the cavity, this simple 
classification is changed in general and we observe ohmic modes and combined 
junction/ohmic modes. The existence or not of the nth resonance is connected 
to the value of the Fourier coefHcient c'^ = cosmra/l. 

These results carry over to the case of an array with many junctions where 
it is possible to choose a voltage where one of the junctions is inactive and 
another is in junction mode. We also use the analysis to understand the eflfect 
of an external magnetic field and the infiuence of having or tt junctions. We 
believe that such a device composed of normal and tt junctions placed at specific 
locations in a microstrip can be used to great advantage for specific applications 
like resonators and frequency mixers. 
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10 Appendix 

10.1 Continuum limit. 

We integrate from a — h to a + h 
The first term yields 

pa+h pa+h 

/ 4>ttdx = / (j)tt{a) {x - a)(j)ttx{a) + 0{\x - a\'^)dx 

J a — h J a — h 



The second one 

r-a+h 



(pxxdx ^ (j)x{a + h) - 4)x{a ~ h) 

a — h 



and the third one 



a+h pa+h 

[sm{(j)) + a(l)t + c(j>tt] dx = / sin [(f>{a) + {x - a)(l)x{a) + Oi{\x - a\'^)) 

—h J a—h 

+a4>tia) + a{x - a)(^-j,t(a) + O^i^x - a^) 
+c4)tt{a) + c{x - a)(t>xtt{a) + Oz{\x - a\'^)dx 

a+h 

sin (0(a) -\- {x — a)(j)x{a) + Oi{\x — a]^)) dx 

-h 

+a4>t{a) + ah(t>xt{.a) + c4>tt{a) + ch(t>xtt{a) + Oi{h^) 
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We now expand the sine term 



sin {(j){a) + {x — a)4>x{a)) 



sm{<j)[a)) cos((a:; — a)(l)x{a)) 
+ sin((a; - a)(l)x{a)) cos(0(a)) 



so that 



/ 

J a 



sin {(f>{a) + (x — a)(l)x{a)) 



a — h 




a+h 



sin((/>(a)) cos((a:: — a)(j)x{a)) + sin((a; — a)(f>x{a)) cos{(f){a))dx 



sin(0(a)) 



2 



sin(ft,02;(i)) — 



Consider the limit of this term when h 



0. We have: 



'°'2h (pxia) 



siTi{h(j)x{a)) — > d, 



so that the equation for the point junction is 

(t)x{a + h) - 4>x{a - h) + da{sin(4>{a)) + a0t(a) + c(f>tt{a)) = j, 
which is consistent with the model H19|l . 



The basis of the method is to discretize the spatial part of the operator and 
keep the temporal part as such. We thereby transform the partial differential 
equation into a system of ordinary differential equations. This method allows 
to increase the precision of the approximation in time and space independently 
and easily. In our case the operator is a distribution so that the natural way to 
give it meaning is to integrate it over a volume. We therefore choose as space 
discretisation the finite volume approximation where the operator is integrated 
over reference volumes. The value of the function is assumed constant in each 
volume. 

As solver for the system of differential equations, we use the Runge-Kutta 
method of order 4-5 introduced by Dormand and Prince implemented as the 
Fortran code D0PRI5 by Hairer and Norsett ^Tj which enables to control the 
local error by varying the time-step. 



0tt - 4ixx + daS{x - a) {sm{(f>) + a(j)t + c(j)tt) = j 



10.2 Method of lines 
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We first transform 12UII into a system of first order partial differential equa- 
tions We write ip{x,t) — (t)t{x,t). 



ip{x,t) = (t>t{x,t) 

Tptix,t) = (f)x^{x,t) - 6{x ~ a){daSiii{4>{x,t)) + a'il!{x,t)) + j 



(67) 



with the boundary conditions : 4'x\l — (f'xl-L — 0- 

For simplicity we will describe the implementation of the finite volume dis- 
cretisation in the case of a single junction. We introduce reference volumes Vk 
whose centers we call a;*;, 1 < A; < nn. The discretisation points are placed such 
that the point Xng+i is at the junction, {x^g^i = a). We thus define Xk and Vk 
using the following identities 



Vk 



with {ng + l)hg = 
Vk 



Xk - —,Xk + — 



_hd hd 
Xk 2 ' ~^ 2 



< k < ng + 1 



ng + 1 < k < nn + 1 



with (nn — ng)hd = I — a. Finally at the junction, k ~ ng + 1 

hn 



Vk^ 



1 ■^ng-\-l 



hd 
2 



nn, ng and nd are respectively the total number of discretisation points, the 
number of points to the left of the junction and the number of points to the 
right. 

For a fixed t, we assume (f){x, t) to be constant on each volume V^, so that 



Xk + -. 



[x, t)dx = h(j){xk,t) , with h — hg ot h = hd 



Integrating over Vk yields: 

1. In the linear part of the PDE : < k < nn + 1 and k ng + 1: 



ip{xk,t) = (l)t{xk,t) 
1pt{Xk t) — 'l'^^k+l,t)-'i<t>{xh,i)+'l>{xk~l,t) 



+ 3 



(68) 



with /i = /ig for < /c < ng -I- 1 or /i = /id for fc > ng -I- 1 . We recognize 
the usual discretisation of the second derivative. 
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.f +d(sin((^)+a(^t)=jl 




Figure 1: Left top pannel, a schematic drawing of a narrow 2D microstrip line 
containing two small Josephson junctions, right top panel section at one of the 
junction. The bottom panel presents the equivalent ID electric circuit. 



2. At the junction: k = ng + 1, we obtain 



S{x — a) {da sin((/)(a;, t)) + a(t)t{x, t)) 



da Sm{4>{x„g+i ,t))+ aipt {Xng+l , t) 

So that the final system is: 



1p{Xng+l,t) = (j)t{Xng+l,t) 

,/, (rr +^ — 4 / <P{x„g + 2,t)-<Pix„g+l,t) _ <P (x „ g + 1 ,t) - (p {x „ g ,t) 

Vt\Xng+l,i) — hg+hd\ hg/2 hd/2 

^ {da sin((/)(x„g+i ,t)) + a(j>t {xng+i ,t)) + j 



hg-\-hd 



(69) 



The ODE system (jfiSI I69|l is then integrated numerically using the D0PRI5 
integrator. 
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Figure 2: plot of (pait) for j = 0.137414 for a junction position a = 3.43, da = 1 
and length I = 10. The left panel shows the ohmic mode regime and the right 
panel the junction mode regime. The numerical solution of 1)20(1 is shown with 
the crosses while the analytical estimates ((33|l and (|35ll are in dashed line. 




Figure 3: Instantaneous voltage (f)t {x, t) as a function of x for different times for 
the ohmic mode: left only oscillation for the junction shown in Fig. [21placed at 
a — 3.43 (left panel) and left and right oscillation for a junction placed dXa — 2 
(right panel). The average voltage is indicated in both figures, V = jl = 1.37 
in the left panel (same as in Fig. and V — ^ {= jl) on the right panel. The 
other parameters are a = 1, and da = I. 
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Figure 4: Instantaneous voltage 4it{x, t) in the junction mode as a function of x 
for 60 values of time equidistributed between and 10. The parameters are the 
same as for Fig. 1. The average voltage V is indicated on the vertical axis and 
the junction's position is indicated on horizontal axis. 
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Figure 5: Two IV characteristic for a single junction placed in the center of the 
microstrip at a = 5 (left panel) and near the right edge at a = 9 (right panel). 
Only the positions of the junction modes are indicated in units of j^. The other 
ticks correspond to the ohmic modes. For each panel two curves are presented, 
the top one is for increasing current while the bottom one is for a decreasing 
current showing a clear hysteresis. The minimum value of the current is 0, the 
maximum is 4 and the stepping is ^. The other parameters are da = a = 1. 
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Figure 6: Plot as a function of time, in each panel, from top to bottom of the 
total energy Et, the energy in the passive region Ep, the power dissipated by the 
quasi-particles Pq and the Josephson energy Ej for a centered junction a = 5 
(bottom panels) and an off-centered junction a = 9 (top panels). The left panels 
show the junction mode n = 4 and the right panels the corresponding ohmic 
mode for the same current. The averages of Et and Ej are shown on the y axis. 
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Figure 7: High voltage case: plot of (j){x, t) as a function of x pour 6 consecutive 
times separated by At = 0.1. The two junctions are located at Oi = 2 and 
02 = 6.2. The other parameters are I = 10, j = 1, and a = 0.5. 
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Figure 8: Three types of ohmic modes obtained for two junctions in a microstrip. 
The positions of the junctions are given in the pictures. We plot </>( {x,t) Vx G 
[0, /] for six consecutive times. The parameters of equation (pn|l are A = 0.5 and 
a = 1. All figures have been obtained with the same current jl — 2.6927937. 
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Figure 9: IV characteristic curves for two junctions placed in ai = 5 and 02 = 
6.2. The critical current densities are di = 0.7, ^2 = 0.3 for the left panel, 
0.5, 0.5 for the middle panel and 0.3, 0.7 for the right panel. We plot on all 
figures jl = VV^ + 1, jl = V, jl = diVV^ + 1 + ^2 V" and jl = d2VV^ + 1+diV. 
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Figure 10: Plots of the phases at the junctions cj){ai,t) and (j){a2,t) as a function 
of time for the IV curves in the two left panels of Fig. |2| The left panels are for 
di = 0.7, d2 — 0.3 and the right panels for di = ^2 = 0.5. The top panels show 
the resonance of order 3 while the bottom panels show the resonance of order 4. 
Note how the junctions switch between ohmic mode and junction mode. The 
analytic expressions for the corresponding junction modes (|35|l are shown in 
dashed line. 
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Figure 11: IV curves for two symmetrically placed junctions ai = 2 et a2 = 
8. The scale on the x axis is in units of j. In the left panel the junctions 
are identical while in the right panel the second one is a 7r-junction. Note 
the resonance on the even harmonics 2^ in the left panel and on the odd 
harmonics (^*'+^)'^ jn the right panel. Wc could have obtained the same change 
of IV curves for the same array of standard junctions by an appropriate choice 
of the magnetic field (see text). 
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